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We use the coupled cluster method in high orders of approximation to make a comprehensive 
study of the ground-state (GS) phase diagram of the spin-i J1-J2-J3 model on a two-dimensional 
honeycomb lattice with antiferromagnetic (AFM) interactions up to third-nearest neighbors. Results 
are presented for the GS energy and the average local on-site magnetization. With the nearest- 
neighbor coupling strength Ji = 1 we find four magnetically ordered phases in the parameter 
window J2, J3 G [0, 1], namely the Neel (N), striped (S), and anti-Neel (aN) collinear AFM phases, 
plus a spiral phase. The aN phase appears as a stable GS phase in the classical version of the model 
only for values J3 < 0. Each of these four ordered phases shares a boundary with a disordered 
quantum paramagnetic (QP) phase, and at several widely separated points on the phase boundaries 
the QP phase has an infinite susceptibility to plaquette valence-bond crystalline order. We identify 
all of the phase boundaries with good precision in the parameter window studied, and we find three 
tricritical quantum critical points therein at: (a) (J^ 1 ,-^ 1 ) ~ (0.51 ± 0.01,0.69 ± 0.01) between 



the N, S, and QP phases; (b) (J 2 C2 , J 3 C 



phases; and (c) (J^ 3 , 



J?) 



(0.65 ± 0.02, 0.55 ± 0.01) between the S, spiral, and QP 



(0.69 ± 0.01, 0.12 ± 0.01) between the spiral, aN, and QP phases. 



PACS numbers: 75.10.Jm, 75.50.Ee, 03.65.Ca 



I. INTRODUCTION 



Phase transitions are common phenomena in nature 
and they have long been the subject of intense theoret- 
ical interest. Since they are driven by thermal and/or 
quantum fluctuations, an understanding of them at the 
microscopic level requires a good many-body description. 
Particular interest in recent years has centered on the so- 
called quantum phase transitions that occur in systems 
in their ground state at zero temperature (T — 0) as 
some parameters in the Hamiltonian describing them are 
varied. 

Spin-lattice systems with competing interactions pro- 
vide a particularly rich field in which to study such quan- 
tum phase transitions. The exchange interactions that 
lead to collective magnetic behavior in spin systems are 
clearly quantum-mechanical in origin. One knows too 
that the interplay between reduced dimensionality, frus- 
tration (due to either competing interactions or the un- 
derlying lattice geometry), and strong quantum fluctu- 
ations generates a huge variety of new states of con- 
densed matter beyond the usual states of quasiclassical 
long-range order (LRO). A particularly rich place to ob- 
serve such exotic ground-state (GS) phases is in situa- 
tions where they originate as a result of quantum fluc- 
tuations within a large set of configurations that are de- 



generate at the classical levelfi^ 

The search for such exotic phases that owe their exis- 
tence purely to quantum effects is nowadays one of the 
primary reasons for the study of frustrated quantum spin- 
lattice systems. The interplay between magnetic frustra- 
tion and quantum fluctuations provides a powerful mech- 
anism for disturbing, destabilizing, or even completely 
destroying magnetic order. The search for a genuine 
quantum spin-liquid (QSL) phased— which has no mag- 
netic order and no LRO or long-range correlations of any 
kind, has itself attracted huge theoretical interest ever 
since its first proposal nearly 40 years ago by Anderson^ 
and its subsequent pursuit by, for example, Shastry and 
Sutherland^, and many others following them. 

The GS wave function of a QSL is clearly a superposi- 
tion of a large number of different configurations. An ex- 
ample could be the resonating valence-bond (RVB) state, 
which is a superposition of many short-range singlet va- 
lence bonds. The RVB state was itself first proposed by 
Fazekas and Anderson^ as the GS wave function for the 
spin-i isotropic Heisenberg antiferromagnet (HAFM) on 
the geometrically frustrated, two-dimensional, (2D) tri- 
angular lattice with only nearest-neighbor (NN) interac- 
tions, all of equal strength, although it is now known 
that this proposal is incorrect and that this system is 
magnetically ordered. 
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In a genuine QSL state no symmetry is broken and 
the quantum fluctuations are required to form a many- 
body singlet state that contains no long-range correla- 
tions with respect to any operator, although there may 
be present some form of topological order. Other ex- 
otic quantum paramagnetic (QP) states, which also have 
no magnetic LRO but which break some spatial sym- 
metry with respect to short-range magnetic correlations, 
can also arise. The various QP valence-bond crystalline 
(VBC) solid phases fall into this category. 

As we have noted, a combination of strong quantum 
fluctuations and strong frustration in a spin system pro- 
vides an ideal scenario for the emergence of such novel 
quantum GS phases as the QSL and other QP phases 
discussed above, which do not possess the magnetic LRO 
that typifies the classical GS phases of the corresponding 
models taken in the limit s — > oo of the spin quantum 
number s of the lattice spins. We know that quantum 
fluctuations tend to be largest for the smallest values of 
s, for lower dimensionality D of the lattice, and for the 
smallest coordination number z of the lattice. Thus, for 
spin-i models the honeycomb lattice plays a special role 
since its coordination number (z = 3) is the lowest pos- 
sible for D = 2. Frustration is easily incorporated by the 
inclusion of competing ncxt-nearcst-ncighbor (NNN) and 
possibly also next-next-nearest-neighbor (NNNN) bonds. 
For these reasons such spin- ^ frustrated Heisenberg mod- 
els on the honeycomb lattice have engendered huge the- 
oretical interest.— ~— 

Interest in the honeycomb lattice has been given fur- 
ther impetus by the discovery of a QSL phase in the ex- 
actly solvable Kitaev model^ in which the spin-i parti- 
cles are sited on a honeycomb lattice. Additional interest 
has also emanated from the recent synthesis of graphene 
monolayers^ and other magnetic materials with a honey- 
comb structure. For example, it is likely that Hubbard- 
like models on the honeycomb lattice may well describe 
many of the physical properties of graphene. In this con- 
text it is particularly interesting to note the clear ev- 
idence of Meng et al£5. that quantum fluctuations are 
strong enough to trigger an insulating QSL phase be- 
tween the non-magnetic metallic phase and the anti- 
ferromagnetic (AFM) Mott insulator for the Hubbard 
model on the honeycomb lattice at moderate values of the 
Coulomb repulsion U. This latter Mott insulator phase 
corresponds in the limit U — > oo to the pure HAFM on 
the bipartite honeycomb lattice, whose GS phase exhibits 
Neel LRO. However, higher-order terms in the t/U ex- 
pansion of the Hubbard model (where t is the Hubbard 
hopping term strength parameter) lead to frustrating ex- 
change couplings in the corres pon ding spin-lattice limit- 
ing model (and see, e.g., Ref. [26[), in which the HAFM 
with NN exchange couplings is the leading term in the 
large- U expansion. 

The unexpected result of Meng et al.,— together with 
other related worker— has excited much interest in un- 
derstanding the physics of frustrated quantum magnets 
on the honeycomb lattice. In particular, a growing con- 



sensus is emerging^ ' n i 14 i 16 i 17 i 19 that frustrated spin-i 
HAFMs on the honeycomb lattice exhibit a frustration- 
induced QP phase. It is interesting to note in this 
context that recent experiments^ on the layered com- 
pound Bi 3 Mn 4 0i2(N0 3 ) (BMNO) at temperatures be- 
low its Curie- Weiss temperature reveal QSL-like behav- 
ior. In BMNO the Mn 4+ ions are situated on the sites 
of (weakly-coupled) honeycomb lattices, although they 
have spin quantum number s = |. The successful sub- 
stitution of the s = | Mn 4+ ions in BMNO by V 4+ ions 
could lead to a corresponding experimental realization of 
a spin-i HAFM on the honeycomb lattice. 

Other realizations of quantum HAFMs which exhibit 
the honeycomb structure include magnetic compounds 
such as InNa 3 Cu 2 SbOg^ and InCu 2 / 3 Vi/ 3 3i 2I In both 
of these materials the Cu 2+ ions in the copper oxide lay- 
ers form a spin-^ HAFM on a (distorted) honeycomb lat- 
tice. Other similar honeycomb materials include the fam- 
ily of compounds BaM 2 (X0 4 ) 2 (M=Co, Ni; X=P, As)^ 
in which the magnetic ions M are disposed in weakly- 
coupled layers where they are situated on the sites of a 
honeycomb lattice. The Co ions have spins s = ^ and 
the Ni ions have spins s = 1 . Recent calculations^ of the 
material /3-CU2V2O7 have demonstrated that its proper- 
ties can also be described in terms of a spin- ^ model on 
an (anisotropic) honeycomb lattice. 

Finally, we note that the very recent prospect of being 
able to realize spin-lattice models with ultracold atoms 
trapped in optical lattices^ is likely to make even more 
data available about the quantum phase transitions in 
the models as the exciting possibility opens up in such 
trapped-atom experiments to tune the strengths of the 
competing magnetic bonds, and hence to drive the sys- 
tem from one phase to another. 

Recently, we have made a series of studies of the frus- 
trated spin-i J1-J2-J3 model on the honeycomb lattice 
using the coupled cluster method (CCM) complemented 
in some cases with the Lanczos exact diagonalization 
(ED) of small lattices. We have studied various regimes 
in the full (Ji, J2, J 3 ) parameter space of NN (Ji) bonds, 
NNN (J 2 ) bonds, and NNNN (J 3 ) bonds. These include 
(a) the AFM model (i.e., with J\ > 0) in the special 
case where the NNN and NNNN bonds are also AFM 
and have equal strength (J 3 = J 2 = kJi > 0)ji£ (b) 
the ferromagnetic (FM) model (i.e., with J\ < 0) with 
frustrating NNN and NNNN bonds, again in the special 
case where they are both AFM and have equal strength 
{J3 = Ji > 0)^ ( c ) the AFM model (i.e., with Ji > 0) 
in the special case where the NNN and NNNN bonds 
are both FM and have equal strength (J 3 — J2 < 0)^ 
and (d) the AFM model (i.e., with J\ > 0) in the spe- 
cial case where we have frustrating NNN bonds only 
(i.e., with J 2 = xJi > 0; J 3 = 0)2± The aim of the 
present work is to extend the investigation of the phase 
diagram of the full spin-i J\-Ji~J?, model on the hon- 
eycomb lattice in the case where all of the NN, NNN, 
and NNNN bonds are AFM, but no further restriction is 
made except that we limit the parameter space window 



3 



to J2/J1, J3/J1 G [0,1]. 

We briefly outline the structure of the rest of the pa- 
per. The model itself is first described in Sec. UH before 
briefly outlining the CCM formalism that we employ as 
our main calculational tool in Sec. IIIII To aid the reader 
we first give an overview of our main results in Sec. IIV1 
focussing on the phase diagram for the model, before we 
give a detailed presentation and discussion of our results 
in Sec. IV] Finally, we conclude in Sec. lVII with a summary 
and comparison of our results with the work of others. 



II. THE HONEYCOMB MODEL 

The spin-i J1-J2-J3 model on the honeycomb lattice, 
or special cases of it (e.g., when J 3 = J 2 or J 3 = 0) 
have been intensively studied by many authors (see, e.g., 
Refs. and references cited therein). The Hamilto- 

nian of the model is 
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where i runs over all lattice sites, and where j, k, and 
I run over all NN sites, all NNN sites, and all NNNN 
sites to i, respectively, counting each bond once and once 
only. Each site i of the lattice carries a spin-s particle 
represented by an SU(2) spin operator Sj = (sf,sf,sf). 
We restrict ourselves here to the case s = 5. The lattice 
and the exchange bonds are illustrated in Fig. [TJ 

Before discussing the s = ^ version of the model that 
is the topic of the present paper it is useful to consider 
first the classical limit (i.e., s — > 00). Thus the Ji~J 2 - 
J3 model on the honeycomb lattice has six classical GS 
phases in the case where J\ > (as considered here) 
and where J 2 and J3 are arbitrary (i.e., can take either 
sign)<££ The six phases comprise three collinear AFM 
phases, the FM phase, plus two different helical phases 
(and sec, e.g., Fig. 2 of Ref. 0). The three AFM phases 
are the Neel (N) phase, the striped (S) phase and the 
anti-Neel (aN) phase as shown in Figs.rjja), (b), and (d) 
respectively. The S, aN, and N states have, correspond- 
ingly, 1,2, and all 3 NN spins to a given spin antiparallel 
to it. Similarly, if we consider the sites of the honeycomb 
lattice as comprising a set of parallel sawtooth (or zigzag) 
chains (in any one of the three equivalent directions) , the 
S state consists of alternating up-spin and down-spin FM 
chains, whereas the aN state consists of AFM chains in 
which NN spins on adjacent chains are parallel. 

Although at T = there exists an infinite family of 
non-coplanar states degenerate in energy with respect to 
each of the S and aN states, both thermal and quantum 
fluctuations^ favor the collinear configurations. When 
J 3 > the spiral state shown in Fig. QJc) is the stable 
classical GS phase in some region of the parameter space. 
This state is characterized by a pitch vector perpendicu- 
lar to one of the three equivalent J\ bond directions, and 
a single spiral angle defined so that as we move along 



the parallel sawtooth chains [drawn in the horizontal di- 
rection in Fig. HJc)] the spin angle increases by n + <j> 
from one site to the next, and with NN spins on adja- 
cent chains antiparallel. The classical GS energy for this 
spiral state is minimized when the pitch angle takes the 
value 4> = cos~ 1 [i(Ji — 2 J 2 )/( J2 — Js)]- The correspond- 
ing minimum value for the GS energy per spin is then 
given as 
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When <f> — > this spiral state simply becomes the 
collinear N state with a corresponding energy per spin 
given by 
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Clearly, the phase transition between this spiral state and 
the N state is of a continuous nature and the correspond- 
ing phase boundary is given by the equation y = -, 
for g < x < \, where y = J3/J1 and x = J 2 /Ji. 

Similarly, when — > ir this spiral state becomes the 
collinear S state with a corresponding GS energy per spin 
given by 



E$ s 2 , 

-JL = -(J 1 - 2 J 2 - 3 J 3 ) . 



(4) 



The classical spiral and S states undergo a continuous 
phase transition along their common phase boundary y = 
I a; + |, for x > \. Furthermore there is a first-order 
phase transition between the collinear N and S states 
along the boundary line x = i, for y > i. These three 
phases (N, S, and spiral) meet at the tricritical point 
(x, y) = (i, |). We- note too that as x — > 00 (for a fixed 
finite value of y) the spiral pitch angle <ft — > |tt. Thus in 
this limiting case the classical model simply becomes two 
HAFMs on weakly connected interpenetrating triangular 
lattices, with the classical triangular-lattice ordering of 
NN spins oriented at an angle |tt to each other on each 
sublattice. 

When y > (and J\ > 0) the above three states are 
the only classical GS phases. When y < the N state 
persists in a region bounded by the same boundary line 
as above, y = ^x — j, for — 1 < x < g, on which it contin- 
uously meets a second spiral state, and by the boundary 
line y = — 1, for x < at which it undergoes a first- 
order transition to the FM state, which itself is the stable 
GS phase in the region x < — \ and y < — 1. Another 
collinear AFM state, the aN state shown in Fig. Hfd), 
with a GS energy per spin given by 
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becomes the stable GS phase in the region x > 
y < \{x - [x 2 + 2{x - \) 2 ] 1/2 }- On the boundary it un- 
dergoes a first-order transition to the spiral state shown 
in Fig. [He). 




FIG. 1: (Color online) The J1-J2-J3 honeycomb model with 
(c) spiral and (d) anti-Neel (aN) states. The spins on lattice 

Finally, for | < x < \ the spiral state shown in 
Fig-QJc) meets a second spiral GS phase on the boundary 
line y = 0, along which there is a first-order transition 
between the two spiral states. This second spiral phase is 
characterized by a pitch vector that is parallel to one of 
the J\ bond directions and by two spiral angles a and /?. 
Within a unit cell the spin directions deviate from those 
in the N state by an angle a, and as one advances from 
one unit cell to the next there is a twist by an angle j3. 
Both of the two pitch angles of this second spiral phase 
smoothly approach the value zero along the above bound- 
ary with the N state, and the value ir along a second 
boundary curve that joins the points (x,y) = (— s,— 1) 
and (|,0), on which it meets the aN state. Both transi- 
tions are continuous in nature. This second spiral phase 
meets the three collinear states N, aN, and FM at the 
tetracritical point (x,y) = (—5, — !•)• 

In this paper we further our study of the spin-i model 
of Eq. ([1]) on the honeycomb lattice, restricting ourselves 
to the case where all of the bonds are antiferromagnetic 
in nature. Thus, henceforth we set J\ = 1 to set the over- 
all energy scale, and we work here within the parameter 
space window J2, J3 £ [0, 1]. 

Recently we have used the CCM to study the spe- 
cial J1-J2 case of this model (and in this parameter 
window with J\ > and J2 = xJ\ > 0) in which 
J3 = Oi^i We found a paramagnetic plaquette valence- 
bond crystalline (PVBC) phase for x ci < x < x C2 , where 
x Cl w 0.207 ± 0.003 and x C2 w 0.385 ± 0.010. We found 
that the transition at x Cl to the N phase appeared to be 
of a continuous deconfined type (although we could not 
exclude a very narrow intermediate phase in the range 
0.21 < x < 0.24), while that at x C2 to the aN phase 
appeared to be of first-order type. As we noted above 
the aN phase exists in the classical version of the J\-J% 
model only at the isolated and highly degenerate critical 
point x = 5. The spiral phases that are present classi- 
cally for all values x > | were found to be absent for all 
x < 1. 

We have also separately used the CCM to study the 
spin-i model of Eq. (JXJ) on the honeycomb lattice in the 
special case where J 3 = J 2 = kJ\ > and Ji > (i.e. 
along the line y = x = k).— We also found a PVBC phase 
in this case for k Ci < K < k C2 > where K Cl ~ 0.47 and k C2 ~ 
0.60. Once again, the evidence favored the transition at 
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Ji = 1; ,h > 0; J 3 > 0, showing the (a) Neel (N), (b) striped (S), 
sites • are represented by the (red) arrows. 

K C1 to the N phase to be of a continuous deconfined type, 
while that at k C2 to the S phase appeared to be first-order 
in nature. 

In order to shed more light on the J1-J2-J3 HAFM 
model (i.e., with Ji > 0) on the honeycomb lattice, we 
now extend its study to map out its entire phase dia- 
gram in the parameter regime x,y G [0,1]. In view of 
the proven success of using the CCM on the above spe- 
cial cases of this model, we continue to utilize it, and in 
Sec. lIIII we accordingly first briefly outline the method as 
it is applied here. 



III. THE CCM FORMALISM 

The CCM is a widely used microscopic many-body 
technique. It is has been demonstrated to be particu- 
larly efficient and very accurate in handling a wide va- 
riety of highly frustrated quantum magnets. Such frus- 
trated systems are notoriously challenging at the theo- 
retical level. Only a very limited number of established 
numerical methods exist for their accurate treatment, 
and other recent and very promising approaches such as 
those based on projected entangled pair states^ have not 
yet been sufficiently widely tested to be able to evaluate 
properly their accuracy and efficacy. 

Of the other well established and widely used tech- 
niques we note that exact diagonalization (ED) tech- 
niques, which involve the finite-size extrapolation of 
numerical exact data for finite-lattice systems, are 
much more challenging for the present honeycomb-lattice 
model than for comparable square-lattice models, to 
which they have been ver y ef ficiently and accurately ap- 
plied (see, e.g., Refs. [39 - [38| ). The reasons include the 
facts that for the honeycomb lattice the unit cell now 
contains two sites, and that there exist relatively fewer 
finite-sized lattices (than in the square-lattice case) that 
are small enough for ED techniques to be used but which 
also contain the full point-group symmetry^ Also, quan- 
tum Monte Carlo (QMC) methods are severely restricted 
in the presence of frustration by the well-known "minus- 
sign problem" . 

By contrast the CCM, when evaluated to high orders 
in one of its systematic approximation hierarchies, as de- 
scribed below, has been proven through a huge variety 
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of applications to provide a powerful tool both to deter- 
mine with good accuracy the positions of quantum crit- 
ical points , 15 i 20 ~ 22 i 38 ~ — and to classify the nature of any 
QP phases in the systemi 15 ' 21 ' 44 ' 48 As noted above, we 
have also used the CCM very successfully in some pre- 
vious applications to honeycomb-lattice models,—" 2 ^— 
and for all these reasons we now employ it again here. 

The CCM is a size-extensive method, in which the limit 
N — > oo, where N is the number of lattice spins, may 
automatically be imposed from the outset. The many- 
body system under study is assumed to have exact ket 
and bra GS energy eigenvectors, and (*B\ respectively, 
which satisfy the corresponding Schrddinger equations, 

H\if>) = E\V) , (^\H = E{^\, (6) 

and which are chosen to have the normalization (\f r |\f r ) 
1 . i.e., (#| = (* 1/^(^1^}. The quantum correlations 
present in the exact ground state are expressed system- 
atically in the CCM with respect to some suitable nor- 
malized model (or reference) state, It is com- 
mon practice to choose simple quasiclassical states as 
CCM reference states, although other choices are cer- 
tainly possible. In this study we choose various classical 
model states as our CCM model states, namely: (a) the 
Neel, (b) the striped, (c) the spiral, and (d) the anti-Neel 
states shown in Fig. [TJ as we discuss below. 

The model state |4>) is required to be a fiducial vector 
in the sense that all possible ket states in the many-body 
Hilbert space can be obtained by acting on it with an 
appropriate linear combination of mutually commuting 
many-body creation operators, Cf , which may be de- 
fined with respect to the model state. The operators 
Cj = (Cj)^, with Cq = 1, thus have the property 
that ($|C+ = = C7|$) = 0; V/ f 0. The CCM 
parametrizations of the exact ket and bra GS wave func- 
tions are given in terms of the usual exponentiated forms, 

|*)=e s |$), <*| = <$|£e- s , (7) 

where the CCM correlation operators S and S arc them- 
selves expressed as generalized multiconfigurational cre- 
ation and destruction operators respectively, 

5 = ^5/C+, Si = l + J2^ c 7 ,VJ^0. (8) 

i i 

Clearly these parametrizations satisfy the normalization 
relations (*|*) = ($1*) = ($|$) = 1- 

The set of correlation coefficients (Si, Si) is now 
determined by requiring the energy expectation value 
H = (\0 , |i/|\I f ) to be a minimum with respect to each of 
the correlation coefficients themselves. This will result in 
the coupled sets of equations [§\C~f e~ s He s '|$) = and 
(®\S(e- s He s - E)C}\<$>) = 0; V/ ^ 0, which we nor- 
mally solve for the correlation coefficients (Si, Si) using 
parallel computing routines once the specific truncation 
scheme is specified, as described further below. 



TABLE I: Number of fundamental configurations, Nf, for the 
spin- i J1-J2-J3 model (Ji = 1) on the honeycomb lattice, 
using the Neel, striped, anti-Neel, and spiral states. 



Method 


Nf 




Neel 


striped 


anti-Neel 


spiral 


LSUB4 


5 


9 


9 


66 


LSUB6 


40 


113 


85 


1080 


LSUB8 


427 


1750 


1101 


18986 


LSUB10 


6237 


28805 


17207 


347287 



In order to treat each lattice site in the spin system on 
an equal basis it is extremely convenient to rotate the lo- 
cal spin-axes on each site in such a way that all the spins 
of each CCM reference state used point along the nega- 
tive z-direction. Such rotations in spin space are obvi- 
ously canonical transformations that have no effect on the 
fundamental SU(2) commutation relations. The spins of 
our system are then represented entirely by these locally 
defined spin coordinate frames. The multispin creation 
operators may be written as linear sums of products of 
the individual spin raising operators = + is v k , i.e., 
Cf = sjj" Sfc ■ • • Sfe • After calculation of the correlation 

coefficients (Si, Si), we can then calculate the GS en- 
ergy using E = ($|e _ He s \<&), and the magnetic order 
parameter, which is defined to be the average local on-site 
magnetization, M = — jj(&\ J2iLi s i with respect to 
the local rotated spin coordinates described above. 

If we include all possible multispin configurations for 
the calculation of the correlation coefficients (Si, Si), 
then the CCM formalism becomes exact. Of course in 
practice one needs to truncate the set, and there are sev- 
eral well-developed and systematically improvable trun- 
cation hierarchies that have been extremely widely tested 
by now. For spin-i systems we usually use the well- 
established localized LSUBm truncation scheme where 
we keep at a given truncation level specified by the trun- 
cation index m only all of those multi-spin configurations 
which may be defined over all possible lattice animals 
(or polyominos) of size m on the lattice. A lattice an- 
imal (or polyomino) of size m is defined as a set of m 
contiguous sites in the usual graph-theoretic sense where 
every site is adjacent (in the nearest-neighbor sense) to 
at least one other site. The method of solving for higher 
orders of LSUBm approximations is well documented in 
Refs. [49ll5l"l ], to which the interested reader is referred 
for further details. 

Tabic U shows the number Nf of fundamental config- 
urations that are inequivalent after all space and point- 
group symmetries of both the Hamiltonian and the model 
state have been taken into account, for each of the Neel, 
striped, spiral, and anti-Neel model states of the spin- 
\ J1-J2-J3 model on the honeycomb lattice. We note 
that the number Nf of such independent spin configu- 
rations taken into account in the CCM correlation op- 
erators S and S increases rapidly with the truncation 
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index to. Clearly the number of independent configura- 
tions is smaller for states such as the Neel state that have 
a higher degree of point-group symmetry, and for which 
we can utilize conservation laws such as sf^ = 0, where 
s T = J2iLi s i 1S the total spin operator referred to the 
global spin coordinates. 

Clearly the spiral state has the largest number Nf, 
for a given level of LSUBto approximation, from among 
our four model states, and for this state we are limited 
to values to < 10 even with the use of massively parallel 
computing to derive and solve the corresponding coupled 
sets of CCM bra- and ket-state equations^ For the spiral 
state we note too that we have the additional computa- 
tional cost that the pitch angle at a given LSUBto level 
must be chosen to minimize the corresponding estimate 
for the GS energy. For specified values of each of the 
exchange parameters J2 and J3 (with J\ = 1) a typical 
computational run for the spiral phase at the LSUB10 
level typically requires about 5 h computing time using 
3000 processors simultaneously. 

Although the CCM works from the outset in the limit 
N — > 00 of an infinite number of spins, and hence the 
need for any finite-size scaling is obviated, we do still 
need to extrapolate the LSUBto data to reach results in 
the exact to —> 00 limit. Although there are no known 
exact extrapolation rules, by now there exists a wealth of 
empirical experience in extrapolating the GS energy, E, 
and the magnetic order parameter (i.e., the average local 
on-site magnetization), M. For the GS energy per spin, 
E/N, a well-established and very accurate extrapolation 



ansatz (see, e.g., Refs. 



IS 



E(m)/N = ao + ftiTO. 2 + a 2 TO 



(9) 



whereas for the magnetic order parameter, M, we use 
different schemes depending on different circumstances, 
specifically on whether the system is highly frustrated or 
not. Thus, for systems with a GS order-disorder transi- 
tion or with a considerable degree of frustration, such as 
is the case for t he pre sent model, we use (see, e.g., Refs. 
[TH3-111.38.42 iimi) 



M(to) = c + CiTO 1/2 + c 2 to" 3/2 . 



(10) 



When we have only three data points to fit to an ex- 
trapolation formula, such as will sometimes occur here, 
specifically for the spiral phase, a two-term extrapolation 
fit can easily be preferable in practice to a three-term fit . 
This is particularly the case when one of the data points 
is either far from the limiting case or when it does not 
represent all of the features of the system as well as the 
remaining, more accurate points. In such cases we some- 
times use the alternative simpler forms, 



and 



E(m)/N = b + b x rn 



M(m) = d + diTO~ 1/2 



(11) 



(12) 



instead of their counterparts in Eqs. ((9j) and (flQ)) . respec- 
tively. 

Finally we note that since the hexagon is an important 
structural element of the honeycomb lattice, it is prefer- 
able for the extrapolations to use only LSUBto data with 
m > 6, wherever possible. However, especially for the 
spiral phase that is particularly costly of computational 
resource, as we explain below, we sometimes need to in- 
clude LSUB4 results in the extrapolations. Under such 
circumstances, however, we always perform a sensitivity 
analysis, by doing some LSUBm runs with higher values 
of m for a few indicative points only, as we discuss in 
more detail in Sec. \V\ 



IV. PREVIEW OF THE PHASE DIAGRAM 

Before discussing our results in detail it is perhaps use- 
ful to summarize our main findings first, and for that 
purpose we show in Fig. [5] the phase diagram for the 
frustrated spin-i model of Eq. (fTJ) on the honeycomb 
lattice, in the case where all the bonds are antifcrro- 
magnetic in nature (i.e., J n > 0, n — 1,2,3). Further- 
more, we set J 1 = 1 and restrict ourselves to the window 
< J m < l,m = 2,3. Henceforth we denote x = J2/J1, 
y = J3/ 'J2. The actual phase boundaries are determined 
from a variety of information that emerges from our CCM 
calculations, as we now describe briefly and with further 
details given in Sec. |V] 

As we have already noted, we have previously stud- 
ied this model for the two special cases with J3 = J2 
in Ref. fTF 



and with J 3 = in Ref. [21|, and the cor- 
responding CCM results from those papers are included 
in Fig. [2] Firstly, along the line J3 = J2 (i.e., when 
y = x = k) we founcU^ that the system has quasiclassical 
AFM Neel (N) order for k < n Cl 0.47, quasiclassical 
AFM striped (S) order for k > n C2 sa 0.60, and a quan- 
tum paramagnetic phase separating the N and S phases 
for k Ci < k < k C2 . By studying the susceptibility of the N 
and S states to hexagonal plaquette valence-bond crystal 
(PVBC) ordering, we found that the most likely scenario 
was that the intervening state had PVBC order over the 



entire range k c 



< K < K r 



The transition at 



between the PVBC and S GS phases was seen to be of 
first-order type, while that at k = n Cl between the N 
and PVBC GS phases appeared to be a continuous one. 
Since the N and PVBC phases break different symmetries 
our results favored the transition point between them at 
K = K Cl to be a deconfincd quantum critical point (QCP). 



The QCPs at y 



and at y 



are clearly 



shown in Fig. [2] with the larger (red) times (x) and the 
larger (green) plus (+) symbols respectively. 

Secondly, in a separate study along the line y = 0, we 
found^i that the system has the quasiclassical N state 
as its GS phase for x < x Cl « 0.21, the quasiclassical 
anti-Neel (aN) state as its GS phase for x > x c . 2 w 0.39, 
and again a quantum paramagnetic (QP) phase separat- 



ing the N and aN phases for x c 



< x < x c 



Similar 
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FIG. 2: (Color online) Phase diagram of the spin-| J1-J2-J3 
model on the honeycomb lattice (with Ji = 1), in the parame- 
ter window J2, J3 G [0, 1]. The five regions correspond to four 
quasiclassical phases with (a) AFM Neel (N) order as shown 
in Fig. []Ja), (b) collinear AFM striped (S) order as shown 
in Fig. [ljb), (c) spiral order as shown in Fig. 0Jc), (d) AFM 
anti-Neel (aN) order as shown in Fig. Did), plus (e) a mag- 
netically disordered, or quantum paramagnetic (QP), phase 
that exhibits plaquette valence-bond crystalline (PVBC) or- 
der on at least part of the boundary region (and see below). 
The first-order phase transition boundary between the N and 
S phases, marked by the (grey) convolution (eight-pointed 
star, ^r) symbols is found from points at which the curves 
for the magnetic order parameter M of the two phases cross; 
the first-order phase transition boundary between the S and 
QP phases, marked by (green) plus (+) symbols, is found 
from points at which M — ► for the S phase; the first-order 
phase transition boundary between the S and spiral phases, 
marked by (cyan) open triangle (A) symbols, is found from 
points at which the curves for the magnetic order parameter 
M of the two phases cross; the phase transition boundary 
between the spiral and QP phases, marked by (orange) open 
circle (O) symbols (of two sizes, see main text in Sec. IV B) . 
is found from points at which M — » for the spiral phase; 
the first-order phase transition boundary between the spiral 
and aN states, marked by (magenta) times (x) symbols, is 
found from points at which the curves for the magnetic order 
parameter M of the two phases cross; the phase transition 
boundary between the aN and QP phases, marked by (blue) 
plus (+) symbols, is found from points at which M — > for 
the aN phase; and the phase transition boundary between the 
N and QP states, marked by (red) times (x) symbols, which 
is probably of continuous (second-order, and possibly of a de- 
confined) nature, is found from points at which M — > for 
the N phase. Points marked by the larger (red) times (x) and 
(green and blue) plus (+) symbols are found to be infinitely 
susceptible to PVBC order, and hence the QP state at these 
points is PVBC in nature. 



CCM calculations of the susceptibility of the N and aN 
phases to PVBC order led again to the conclusion that 
the transition between the PVBC and aN phases was of 
first-order type, while the likely scenario for the transi- 
tion between the N and PVBC phases is again that it 
is of the continuous dcconfined type. Nevertheless, due 



to the difficulty in determining the lower critical value of 
x at which PVBC order is established as accurately as 
we determined the value x = x Cl at which Neel order is 
destabilized, we could not exclude a second scenario in 
which the transition between the N and PVBC phases 
proceeds via an intervening phase (possibly even of an 
exotic spin-liquid variety) in the very narrow window 
0.21 < x < 0.24. Again, the QCPs at (x,y) = (x Cl ,0) 
and (x C2 , 0) are clearly shown in Fig. [2] with the larger 
(red) times (x) and the (blue) plus (+) symbols respec- 
tively 

We also showed previously^ by a comparison of the 
GS energies of the spiral and aN phases calculated sep- 
arately with the CCM, that the spiral phases that are 
present classically (i.e., for the case where the spin quan- 
tum number s — > 00) in the case y = for all values x > g 
are absent for all values x < 1. The actual phase bound- 
ary between the spiral and aN phases shown in Fig. [2] is 
now calculated in the present paper, as described below. 

Based on our previous findings for the GS phases of 
these two special cases when (a) J3 = J2 and (b) J3 = 0. 
we have now performed a series of CCM calculations 
based on the N, S, aN, and spiral states as model states, 
for a variety of cuts in the phase diagram at both con- 
stant values of J3 and constant values of J2. For example, 
the phase boundary between the N and the S phases is 
obtained, as explained more fully in Sec. IV Al from our 
extrapolated (m — > 00) LSUBm results for the order pa- 
rameter M (namely, the average onsitc magnetization) 
of the two phases, for a variety of constant J3 cuts. We 
find that for values of y = J3/J1 > 0.69 the two magneti- 
zation curves meet at a (positive) nonzero value, indica- 
tive of a direct first-order transition between the states. 
These points are shown in Fig. [2] by the (grey) convolu- 
tion (eight-pointed star, symbols. 

For the value y w 0.69 the two curves become zero at 
precisely the same point, x ~ 0.51. Conversely, when 
y < 0.69, the order parameters of the N and the S 
phases both become zero at respective critical values of 
x before the curves cross (when solutions exist for both 
phases), indicating the emergence of a new phase separat- 
ing them. The corresponding points where the magnetic 
order parameters for Neel order and striped order vanish 
are shown in Fig. [2] by (red) times ( x ) and (green) plus 
(+) symbols respectively. By continuity with our earlier 
results^ along the line y — x, we tentatively identify 
the intervening phase as the PVBC state. The tricrit- 
ical QCP between the N, S, and PVBC phases is thus 
identified as being at (x,y) « (0.51,0.69). We also note 
that for values of y < 0.55 no solution for the S phase 
exists with M > for any value of x, giving preliminary 
indications of a new phase boundary between the S state 
and another phase that we identify as a spiral phase. 

By comparing the order parameters for the S and spi- 
ral phases at various constant J2 cuts we find that for 
values of x > 0.66 the two curves meet at a (positive) 
nonzero value, once again indicative of a direct transition 
between the states. These points are shown in Fig. fj] as 
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(cyan) open triangle (A) symbols. For the value x « 0.66 
the two curves become zero at the same point y « 0.55. 
Then, for values x < 0.66 the order parameters of the S 
and spiral phases both become zero at respective criti- 
cal values of y before the curves cross. Once again this 
indicates a phase separating the S and spiral phases for 
values of x < 0.66 (down to a lower value of x w 0.635 be- 
low which the spiral phase ceases to exist for any value of 
y), which we similarly identify tentatively as the PVBC 
phase. 

In that very narrow window 0.635 < x < 0.66, which 
is almost certainly an artifice of our approximations, we 
denote in Fig. [5] the points where the magnetic order 
parameter vanishes (M — > 0) for the striped and spiral 
states by (cyan) open triangle (A) and (orange) open 
circle (O) symbols respectively. We argue in Sec. IV CI 
that these results arc consistent with the existence of a 
second tricritical QCP at (x, y) « (0.65, 0.55) between 
the S, spiral, and (tentatively) PVBC phases. The re- 
mainder of the phase boundary between the spiral and 
PVBC states is similarly identified by the vanishing of 
the magnetic order parameter of the spiral phase, and 
these points are again shown in Fig. [2] as (orange) open 
circle (O) symbols. 

Finally by comparing the energies of the aN and spi- 
ral phases we find that for all values of the parameter 
J2 < 1 where the spiral phase exists, the aN phase ac- 
tually has a lower energy for values of the parameter 
J3 below a certain critical value, which itself depends 
on J-2- Similarly, by comparing the order parameters of 
these two phases at various constant Ji cuts, we find 
that for values of x > 0.69 the two curves meet at a 
(positive) nonzero value, indicative once more of a di- 
rect phase transition between the aN and spiral phases. 
These points are shown in the phase diagram of Fig. [5] 
by (magenta) times (x) symbols. 

For the value x » 0.69 the two curves become zero at 
the same point y w 0.12. Conversely, for values x < 0.69 
the order parameters of the aN and spiral phases both 
become zero at respective critical values of y before the 
curves cross. This is again indicative of a phase sepa- 
rating the aN and spiral phases for values of x < 0.69 
(down to the lower value of x sa 0.635 below which the 
spiral phase ceases to exist for any value of y), as noted 
above. This intermediate phase is again tentatively iden- 
tified as having PVBC order. In this way we identify 
a third tricritical QCP at (x,y) «a (0.69,0.12) between 
the spiral, aN and (tentatively) PVBC phases. Remain- 
ing points on the phase boundary between the aN and 
PVBC phases are then identified as the points where the 
magnetic order parameter of the aN phase vanishes, and 
these are shown by (blue) plus (+) symbols on the the 
phase diagram of Fig. [51 

In Sec. [V] we now describe in more detail how the var- 
ious points in the phase diagram of Fig. [5] arc obtained. 
We also discuss the properties of the various phases that 
we have examined. 



V. RESULTS 

In this section, we present and discuss our CCM results 
for the spin-i J1-J2-J3 HAFM on the honeycomb lattice, 
with all of the bond strengths positive (i.e., antifcrromag- 
netic in nature). To set the overall energy scale we put 
J\ = 1, and we investigate the parameter space window 
^2,^3 £ [0,1]- We use each of the Neel (N), collincar 
striped (S), spiral, and anti-Neel (aN) states shown re- 
spectively in Figs, [lja)-(d) as CCM model states. 



A. Neel versus striped phases 

Figurcs[3Ua) and[3](b) show the extrapolated (to — > 00) 
CCM LSUBto values for the GS energy per spin for the 
Neel (N) and striped (S) states as functions of J2 for 
various fixed values of J 3 in the range 0.5 < J3 < 1.0. 
The extrapolations have been performed using Eq. ([9]) 
and the calculated LSUBto results with to = {6,8, 10}. 
We observe that the energy curves cross for all values 
of the parameter J3 > 0.68, but for values J3 < 0.68 the 
curves do not cross. This gives us a first indication of the 
emergence of an intermediate phase between the N and S 
states, over a finite range of values of the J2 parameter, 
below some critical value of the J3 parameter. 

We note that the extrapolations become more difficult 
in the vicinity of this critical point, and consequently the 
actual values of J2 at which the curves cross for fixed 
values of J3 near the critical value are more uncertain 
than those at larger values. Furthermore, at the actual 
energy crossing points very near the critical point the 
corresponding values of the magnetic order parameter 
(i.e., the average onsite magnetization) M for one or both 
states becomes negative and hence unphysical. Indeed, 
for the S state, M < for the entire J3 = 0.5 curve, 
which is why we have not shown it in Fig. Eta). 

In order to obtain more accurate values of the critical 
point we also show in Figs.@|a) and0|b) the curves for 
the extrapolated order parameters M of the N and S 
states, corresponding to values of J3 shown in Figs. [3^a) 
andEtb) for the GS energy per spin, E/N. The LSUBoo 
curves shown use the LSUBto results with to = {6, 8, 10}, 
together with the extrapolation scheme of Eq. (p~0|) . which 
is appropriate for this highly frustrated regime. 

We observe again that the curves intersect for values 
•h <; 0.69, and that the corresponding values of (J 2 , J 3) 
are our best estimate for the phase boundary between 
the N and S states, as shown on the phase diagram of 
Fig. [3 by the points denoted with (grey) convolution 
(eight-pointed star, *) symbols. For values J3 < 0.69 
the extrapolated order parameters of both the N and S 
phases become zero before the curves intersect, revealing 
the presence of an intermediate phase in that regime. The 
corresponding points in the case J3 < 0.69 where M — > 
for the N and S phases are shown in the phase diagram 
of Fig.[2]by (red) times (x) and (green) plus (+) symbols 
respectively. Our best value for the corresponding tricrit- 
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FIG. 3: (Color online) Extrapolated CCM LSUBoo results for the GS energy per spin, E/N , as a function of J 2 , for various 
fixed values of J3 in the range 0.5 < J3 < 1.0, for the Neel and the striped states of the spin-^ J1-J2-J3 model on the 
honeycomb lattice (with Ji = 1). The extrapolated LSUBm (m — > 00) results are based on the extrapolation scheme of Eq. @ 
and the calculated results with m = {6, 8, 10}. 
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FIG. 4: (Color online) Extrapolated CCM LSUBoo results for the GS magnetic order parameter, M, as a function of J2, for 
various fixed values of J3 in the range 0.5 < J3 < 1.0, for the Neel and the striped states of the spin-i J1-J2-J3 model on 
the honeycomb lattice (with Ji = 1). The extrapolated LSUBm (m — )■ 00) results are based on the extrapolation scheme of 
Eq. (110[) and the calculated results with m = {6, 8, 10}. 



ical QCP comes from the data shown in Fig.^Jb), where 
it is seen to be at (J 2 C \ J 3 Cl ) = (0.51 ± 0.01, 0.69 ± 0.01), 
and where the error bars are estimates from a sensitivity 
analysis of the LSUBm extrapolation scheme. 

We note that the extrapolated order parameter M be- 
comes everywhere negative (i.e., for all values of J 2) for 
the S state for all values of J3 < 0.55, as may be seen 
from data similar to those shown in Fig. E^a). This is 
a clear first indication that the S state becomes unsta- 
ble as the GS phase in this regime. From a comparison 
with the corresponding classical model (i.e., in the limit 
s — > 00) discussed in Sec. [TT1 we might expect the S state 
to yield to the spiral state, at least for sufficiently large 
values of J 2 in the present s = \ case. We investigate 
this further in Sec. IV Cl below. It is clearly also expected 
that the Neel (N) phase will not survive for large enough 
frustrating values of J2 > 0, and again from a compar- 



ison with the classical model we expect that the spiral 
phase might exist in that case too. Hence, we first make 
a comparison in Sec. IVBl of the N and spiral phases. 



B. Neel versus spiral phases 

We start by analyzing the GS energy per spin, E/N , 
for the spiral state as a function of the spiral pitch angle, 
<\>. In our CCM calculations we choose the angle <j>, for 
each point in the phase diagram where the spiral state ex- 
ists, as the one that minimizes the energy estimate there. 
Clearly the minimizing angle in general also depends on 
the particular LSUBm approximation being used with 
the spiral state as CCM model state. For example, we 
show in Fig. [5][a) the GS energy per spin, E/N, as a 
function of pitch angle (f> in the LSUB6 approximation, 
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FIG. 5: (Color online) (a) LSUB6 results for the GS energy per spin of the spin-i J1-J2-J3 model on the honeycomb lattice 
(with Ji = 1), using the spiral state as CCM model state, as a function of the spiral pitch angle <j), for some illustrative values 
of J2 in the range < J2 < 1 and for a fixed value of J3 = 0.4. (b) The pitch angle cj> = ^lsub™ that minimizes the energy 
E^svBm(4>) of the spiral state of the spin-i mo del on the honeycomb lattice (with Ji = 1). The CCM LSUBm results 

with m — 8 are shown as functions of J2 for several fixed values of J3 in the range 0.2 < J3 < 0.7. The Neel results with cj> — 
are also shown artificially offset above the actual <j> — axis for the sake of clarity. 



for various illustrative values of J2 at a fixed value of 
J 3 = 0.4 (and Ji = 1). 

We note first from Fig.O^a) that for various fixed values 
of J2 CCM solutions at a given LSUBm level of approx- 
imation exist only for certain ranges of the spiral angle 
<f>. For example, for J 2 = (and J 3 = 0.4), the CCM 
LSUB6 solutions based on the spiral state exist only for 
< 4> < 0.127T. In this case, where the Neel state (i.e., 
where = 0) is the stable GS phase that minimizes the 
energy, if we try to force the system too far away from 
Neel collinearity the CCM equations themselves become 
unstable in the sense that they no longer have a real solu- 
tion. We note too that as Ji is increased slowly (at fixed 
J3 = 0.4), the minimum in the energy curve at (f> = 
becomes shallower, so that by the time J2 = 0.4 it has 
almost disappeared. This is a first indication of the im- 
minent instability of the Neel state as the GS phase if J2 
is increased slightly more. 

Similarly, for J 2 = 1 (and J 3 = 0.4), the CCM 
LSUB6 solutions based on the spiral state exist only for 
0.43 < 4>/n < 1. In this case, a spiral state (i.e., with a 
value (j> 7^ 0, 7r) is the stable GS phase that minimizes the 
energy, and if we now try to force the system too close 
to the Neel regime, the CCM solution collapses. We also 
observe from Fig.[SJa) that for the smaller value J2 = 0.8 
(and J3 = 0.4), while the energy curve still shows a global 
minimum for a noncollinear spiral phase, it has now also 
developed a secondary minimum at a value <f> = n (i.e., 
that of the collinear striped state), which indicates the 
proximity of the phase boundary between the spiral and 
striped states, as we examine more fully in Sec. IVCI bc- 
low. 

Conversely, as J3 is increased further (for fixed J2 ) , the 
spiral minimum becomes more pronounced, and as J 2 — > 
00 the pitch angle <j> — > |7T. This is as expected, since 



in this limit the model becomes two weakly connected 
HAFMs on interpenetrating triangular lattices, with the 
classical ordering of NN spins oriented at angles |7r with 
respect to one another on each sublattice. 

From data such as that shown in Fig. G3[a) we can cal- 
culate in a given LSUB?n approximation based on the spi- 
ral state as the CCM model state, the angle <j> = 4>hSVBm 
that minimizes the energy, -ELSUBm(</>) for given values 
of the exchange coupling strengths. For example, in 
Fig. E](b) we show the angle 4> = ^lsubs from the LSUB8 
approximation, as a function of the parameter J 2 for sev- 
eral fixed values of the parameter J3 . There is clear pre- 
liminary evidence that for values of J3 below some upper 
critical value there is no stable spiral solution for any 
value of </> 7^ over a certain range of the parameter J2, 
which itself depends on J3. 

Thus, we are led to expect a second tricritical QCP 
in the (J2, J3) plane at (^ 2 2 ,^3 2 ), with J3 2 < J3 1 , such 
that: (a) for values J3 > J3 1 the N and S states meet at a 
common phase boundary discussed in Sec. IV Al above. (b) 
for values J3 1 > J3 > J3 2 there is an intermediate phase 
between the N and S states, and (c) for values J3 < J' 3 2 
there is an intermediate phase between the N state and 
the spiral state with cf> ^ tt. Thus, the QCP at (Jj 2 , J3 2 ) 
is a tricritical point between the N, S, and intermediate 
phases. From the results discussed in Sec. IV Al we now 
expect that J3 2 w 0.55, and we discuss this further below. 

Figure [6] shows our extrapolated CCM LSUBoo results 
for the GS energy per spin, E/N, as a function of J2, for 
various fixed values of J3 in the range 0.2 < J3 < 0.6, 
for the Neel and the spiral states. Once again, the figure 
clearly illustrates the existence of an intermediate phase 
between the Neel and the spiral phases (including the 
striped state as a special case of the latter) for values 

J 3 < J 3 C1 - 
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FIG. 6: (Color online) Extrapolated CCM LSUBoo results 
for the GS energy per spin, E/N, as a function of J2, for 
various fixed values of J3 in the range 0.2 < J3 < 0.6, for the 
Neel and the spiral states of the spin-i J1-J2-J3 model on the 
honeycomb lattice (with Ji = 1). The extrapolated LSUBm 
(m — > 00) results are based on the extrapolation scheme of 
Eq. ^ and the calculated results with m = {6, 8, 10} for the 
Neel state, and with m — {4, 6, 8} for the spiral state. For 
the spiral state the results use the pitch angle cj> = ^LSUBm 
that minimizes the energy E = i?LSUBm (</>)■ 



FIG. 7: (Color online) Extrapolated CCM LSUBoo results 
for the GS magnetic order parameter, M, as a function of J2, 
for various fixed values of J3 in the range 0.2 < J3 < 0.6, 
for the Neel and the spiral states of the spin-i j 1 -j 2 -j i 
model on the honeycomb lattice (with Ji = 1). The extrap- 
olated LSUBm (m —¥ 00) results are based on the extrap- 
olation scheme of Eq. ()12|) and the calculated results with 
m = {6, 8, 10} for the Neel state, and with m = {4, 6, 8} 
for the spiral state. For the spiral state the results use 
the pitch angle <f> = 0LSUBm that minimizes the energy 
E = Blsub™ (</>)■ 



On a technical point we remark that for the spiral state 
the extrapolations are calculated using the LSUBm cal- 
culated results with m = {4,6,8}, rather than with the 
set m = {6, 8, 10} used for the Neel state. This is partly 
due to the very high number, Nf — 347287, of configura- 
tions needed for the spiral state at the LSUB10 level of 
approximation, compared with the corresponding much 
smaller number, Nf = 6237, for the Neel state, as seen 
from Table HI This difference is compounded by the fact 
that for the spiral state we also need to do LSUBm runs 
for each point in the phase space as a function of the pitch 
angle 0, in order to determine the angle <p = ^LSUBm 
that minimizes the corresponding estimate for the en- 
ergy, E = -ELsuBm^)- This makes LSUBm calculations 
for the spiral state with m > 10 particularly demanding 
of computational resources. 

Nevertheless, we did perform LSUB10 calculations for 
the spiral state for the special case J3 = in our pre- 
vious study of the J1-J2 model^ where we performed 
separate extrapolations using the LSUBm results with 
m = {6,8,10} and m = {4,6,8}. We found that both 
extrapolations were in very good agreement with one an- 
other, and hence now feel confident that the spiral-state 
extrapolations for the full J1-J2-J3 model considered 
here with the limited set m = {4, 6, 8} will be equally 
robust, since it is now prohibitively expensive of compu- 
tational resource to perform LSUB10 calculations for the 
spiral state over the whole region of phase space where 
it is the stable GS phase. 

We note that, as is usually the case, the CCM LSUBm 
results for finite m values for a given phase extend beyond 



the actual physical LSUBoo boundary for that phase. 
Thus, the energy curves shown in Fig. [6] for fixed values 
of J3 terminate at certain values of J2, which are deter- 
mined by the termination points of the highest LSUB?7i 
approximations used in the extrapolations, beyond which 
no real solution exists for the corresponding coupled 
CCM equations. We note from Fig. [5] that the maxima 
in the energy curves occur close to these LSUBm termi- 
nation points for the largest m values employed, which in 
turn lie close to the physical (LSUBoo) phase transition 
points. 

It has been suggested^ that such an energy maximum 
approximately coincides with the avoided level crossing 
to a different phase, in which case it could be taken as 
an approximation for the phase transition point of cither 
(the Neel or spiral) state to the intermediate (as yet un- 
known) state. However, we do not use this criterion here, 
since our results for the magnetic order parameter give 
us much more accurate estimates, as we now discuss. 

Thus, we show in Fig. [7] our extrapolated CCM 
LSUBoo results for the GS magnetic order parameters, 
M, for both the Neel and the spiral states, as functions of 
J2, for the same fixed values of J3 shown in Fig. [6] for the 
GS energy. The extrapolated LSUBjti (m — > 00) results 
for the Neel state are based on the extrapolation scheme 
of Eq. (fit)]) with m = {6, 8, 10}, whereas the extrapolated 
results for the spiral state are based on the extrapolation 
scheme of Eq. (jT2j) with m = {4, 6, 8}. 

As we indicated above, LSUB10 calculations are pro- 
hibitively expensive for the spiral state, and hence we 



12 



need for extrapolation purposes to include the LSUB4 
results. When this point is included, and the data set 
to = {4, 6, 8} is thus employed, it is clearly preferable 
to use the extrapolation scheme of Eq. ([12")) rather than 
that of Eq. (jTDJ) , so as not to give the to = 4 result too 
much weight. However, in order to check our results, we 
have performed LSUB10 calculations for the two values 
Js = 0.2, 0.4. For these two values we have also made 
extrapolations using Eq. (|T2|) with to = {6,8,10}, and 
these are indicated in Fig. [2] by the larger (orange) open 
circles. We find, very gratifyingly, that the two extrap- 
olations agree very well with one another at both values 
J 3 = 0.2, 0.4, which gives credence to our results using 
the data set to = {4, 6, 8} elsewhere for the spiral state. 

We note too that in our earlier study2i of this model 
with J 3 = (and J\ = 1), we also computed LSUB12 
results for the Neel state and found that the Neel order 
vanished at a value J2 ~ 0.207±0.003 when we performed 
extrapolations including the to = 12 point. It is this 
point that is shown in Fig. [2] by the larger (red) times 
(x) symbol, although the value obtained with the more 
limited data set m = {6, 8, 10} used for the results in 
Fig. [7J is in good agreement with it. Similarly, in our 
earlier study of the model along the J3 = J2 line^ we also 
used LSUBto results with to = {6, 8, 10, 12} to perform 
the extrapolations, and found that in this case Neel order 
vanished at a value J2 w 0.466 ± 0.005, and this value is 
also shown in Fig. [2] by a larger (red) times ( x ) symbol. 

From curves such as those shown in Fig. [7] we use the 
points where the extrapolated values of the order param- 
eter M vanish, for various fixed values of J3, to plot the 
phase boundaries of the Neel and spiral phases denoted 
in Fig. [5] by (red) times ( x ) and (orange) open circle 
symbols respectively. As expected from our previous dis- 
cussion in Sec. IV A) a Neel-ordered phase exists for all 
values of J 3 up to some critical value of J2 which marks 
its phase boundary. For values J3 < J3 1 « 0.69±0.01 this 
phase borders a quantum paramagnetic phase, whereas 
for J3 > J3 1 it borders the striped state at a first-order 
phase transition boundary. We also see from curves such 
as those shown in Fig. [7] clear evidence for an interven- 
ing phase between the Neel and spiral phases (with pitch 
angle <j> 7^ 7r) everywhere that the spiral phase exists. 
Instead the spiral phase meets the striped phase along 
a common boundary (on which <f> = 71") for all values 
J2 > J| 2 > 0.65. There is thus a second tricritical point 
at (J2 2 , J3 2 ), as we discuss more fully below in Sec. lVC| 
at which the striped, spiral and quantum paramagnetic 
phases meet. 



C. Striped versus spiral phases 

We first recall that classically (i.e., when s — > 00) we 
have for this model (with Ji = 1) that for a fixed value 
of J2 > I the GS phase is the striped phase for J3 > 
5J2 + \ and the spiral phase for J 3 < | J 2 + |. There is 
a continuous phase transition between the two classical 



states along the boundary line, J 3 = \ J2 + \ , J 2 > \, 
on which the spiral pitch angle <p = n - Our results for 
the present s = | model, as we shall see below, indicate 
that quantum fluctuations tend to stabilize the collinear 
order of the striped state to lower values of J3, for fixed 
J2, than the classical limit. Furthermore, as we shall see, 
the quantum fluctuations also seem to turn the classical 
second-order transition into a quantum first-order one. 

Thus, we show in Fig. [SJa) the angle, <f> = </>lsubs that 
minimizes the energy, -Elsubs {<ft) 1 as a function of J3, us- 
ing the spiral state as our CCM model state, for various 
fixed values of J2 . Very similar curves are found for other 
LSUBto approximations. We observe that, unlike in the 
classical case, where 4> — >• tt continuously at the critical 
value, there is now a discontinuous jump on the phase 
boundary. Its origin lies in the double-minimum struc- 
ture of the corresponding energy curves (for fixed values 
of j r 2 and J3) as functions of the pitch angle 0, compa- 
rable to that shown in Fig. [SJa) for the case J2 = 0.8, 
J 3 = 0.4. 

Clearly, if we consider the angle (j> itself to be an order 
parameter (such that <f> = ~k for striped order and 4> 7^ 
0, 7r for spiral order) the typical scenario for a first-order 
transition, as now seen here, is the emergence of such a 
two-minimum structure for E/N as a function of <j) for 
fixed coupling parameter strengths, one at a value <j) 7^ tt 
and the other precisely at <f> — n. For a fixed value of 
J2 we find, at a given LSUBto level of approximation, 
that when J3 is above a certain critical value the global 
minimum in the E = E(4>) curve is at = n, whereas 
when J3 is below this value the global minimum is at the 
other minimum, <p 7^ it. 

Figure [8jb) shows the corresponding extrapolated 
CCM LSUBoo results for the GS energy per spin, E/N, 
of the spiral and striped states, as functions of the param- 
eter J3, for the same fixed values of J2 shown in Fig.[8ja). 
The first-order transition between the spiral and striped 
states can clearly be seen to occur close to, but not pre- 
cisely at, the corresponding maxima in the energy curves. 

As before, the actual phase boundary is most clearly 
seen from our similarly extrapolated CCM LSUBoo re- 
sults for the magnetic order parameter, M, which are 
shown in Fig. [5J We note that for all values of J2 > 0.66 
there is a clear and sharp minimum in the magnetic order 
parameter at the phase transition point in the parame- 
ter J3 where the striped and spiral phases meet. These 
points arc indicated by the (cyan) open triangle (A) sym- 
bols in the phase diagram shown in Fig. [2] 

At the value J2 ~ 0.66 the two curves meet at M = 0. 
We note that for this value of J2 the magnetic order 
parameter M for the spiral state is very small (and pos- 
itive) for all values of J3, and that as J2 is decreased 
further the spiral state rapidly disappears altogether for 
■h ;$ 0.635. In the very narrow regime 0.635 < J2 < 0.66, 
we see from Fig. [9] that there appears to be an intru- 
sion of the intermediate (quantum paramagnetic) phase, 
as shown in Fig. [5] by the appearance of both (cyan) 
open triangle (A) and (orange) open circle (Q) sym- 
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FIG. 8: (a) The angle <f> = </>LSUBm that minimizes the energy i?LSUBm(0) of the spin-i J1-J2-J3 model on the honeycomb 
lattice (with J\ = 1). The CCM LSUBm results with m = 8 are shown as functions of J3 for several fixed values of J2 in the 
range 0.7 < J3 < 1.0. Note that <j> — n corresponds to the striped state, (b) Extrapolated CCM LSUBoo results for the GS 
energy per spin, E/N, as a function of J3, for various fixed values of J2 in the range 0.7 < J3 < 1.0, for the spiral and the 
striped states of the spin-i j 1 -j 2 -j 3 m odel on the honeycomb lattice (with Ji = 1). The extrapolated LSUBm (m — > 00) 
results are based on the extrapolation scheme of Eq. (|9j) and the calculated results with m = {4, 6, 8}. For the spiral state the 
results use the pitch angle <j> — (j>hSVBm that minimizes the energy E = -ELsuBm(0)- 
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FIG. 9: (Color online) Extrapolated CCM LSUBoo results 
for the GS magnetic order parameter, M, as a function of J3, 
for various fixed values of J2 in the range 0.64 < J2 < 1.0, 
for the spiral and the striped states of the spin-i j 1 -j 2 -j i 
model on the honeycomb lattice (with Ji = 1). The extrap- 
olated LSUBm (m — > 00) results are based on the extrap- 
olation scheme of Eq. (|12p and the calculated results with 
m = {4, 6, 8}. For the spiral state the results use the pitch an- 
gle <j> = <?!>LSUBm that minimizes the energy E = S LSU Bm (</>)■ 



bols at the striped-spiral phase boundary at the two 
values J2 = 0.64, 0.65. It seems almost sure, however, 
that this effect arises from our extrapolations, and is 
an indication of the (small) errors inherent in them. 
Our best estimate from the results shown in Fig. |H] is 
thus that the second tricritical QCP, where the spiral, 
striped and quantum paramagnetic phases meet, occurs 
at (J 2 C2 , J3 2 ) = (0.65 ±0.02, 0.55 ±0.01). 

We also note from Fig. [S] that for values 0.635 < J2 < 
0.77 and J3 > the magnetic order parameter M of the 



striped state becomes zero at a lower critical value of J3. 
These lower values in each case are shown in Fig. [2] by the 
same (orange) open circle (O) symbols as we discussed 
previously in Sec. IV Bl We note that for the special case 
J3 = that we investigated earlier^ the spiral state is 
actually unstable, since the anti-Neel state was seen to 
have lower energy for all values of J2 in the range inves- 
tigated, namely J% < 1, where solutions for the spiral 
state could be found. From continuity, we expect that 
the anti-Neel state should remain the stable GS phase 
for small enough values of J3 below some critical value 
for each fixed value of J2, above which value the spiral 
phase then becomes the stable GS phase. Thus we are led 
to expect that there might exist a third tricritical QCP 
at (J2 3 , J3 3 ) between the spiral, quantum paramagnetic, 
and anti-Neel GS phases. We examine this further in 
Sec. [VD] below. 



D. Spiral versus anti-Neel phases 

In Fig. [10] we show the extrapolated CCM LSUBoo 
results for the GS energy per spin, E/N, of both the 
spiral and anti-Neel states, as functions of the parameter 
J3, for various fixed values of the parameter J2 in the 
range 0.7 < J2 < 1.0. Although the energy differences 
are small for each fixed value of J2, the results at each 
LSUBm level, as well as the extrapolated results, clearly 
show an energy crossing point. These energy crossing 
points are thus our first estimates of the phase boundary 
points between the spiral and anti-Neel states. 

On a technical point, we have noted previously that 
CCM LSUBm calculations for the spiral state are com- 
putationally very expensive for values of the truncation 
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FIG. 10: (Color online) Extrapolated CCM LSUBoo results 
for the GS energy per spin, E/N, as a function of J3, for vari- 
ous fixed values of J2 in the range 0.7 < J2 < 1.0, for the anti- 
Neel and the spiral states of the spin-i J1-J2-J3 model on the 
honeycomb lattice (with Ji = 1). The extrapolated LSUBm 
(m — > 00) results are based on the extrapolation scheme of 
Eq. (|11[) . and the calculated results with m = {4,6,8} in 
both cases. For the spiral state the results use the pitch angle 
4> = 0LSUBm that minimizes the energy E = i?LSUBm(^)- We 
note that in all cases curves without symbols attached refer 
to the anti-Neel state, whereas the corresponding curves with 
symbols refer to the spiral state. 



index to > 10. Thus, we are for the most part restricted 
to the data set to = {4, 6, 8} for the spiral state, although 
we have performed a very few calculations for a few se- 
lect values in the parameter space with to = 10. With 
only three data points to fit to an extrapolation formula, 
a two-term extrapolation fit (such as those in Eqs. (TIT]) 
and (Ti"2]) , for example) can often be preferable in practice 
to a three-term fit (such as their counterparts in Eqs. Q 
and (jTUJ) , respectively) . This is particularly the case when 
one of the data points is either far from the limiting case 
or when it does not represent all of the features of the 
system as well as the remaining, more accurate points, 
as is the case here for the to = 4 points. 

Thus, since the energy differences of the spiral and aN 
states are relatively small, we have found it preferable 
to use the extrapolation scheme of Eq. (TIT]) in this case, 
and to employ the same data set with m = {4,6,8} for 
both (aN and spiral) phases, even though results with 
to = 10 are more readily available for the aN state. 
We have, however, demonstrated that the results so ob- 
tained arc robust and reliable, by making further checks 
in some limited test cases using the extrapolation scheme 
of Eq. © fitted to the results to = {4,6,8,10} or the 
extrapolation scheme of Eq. (jlll) fitted to the results 
m = {6, 8, 10}, for example. 

We note that real CCM LSUBm solutions based on the 
aN state as model state cease to exist, for a fixed value 
of J2 , above some termination value in the parameter J3 
that itself depends on the truncation index to, just as 
we have indicated above for other phases. These LSUB8 
terminations are what cause our extrapolations for the 
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FIG. 11: (Color online) Extrapolated CCM LSUBcxj results 
for the GS magnetic order parameter, M, as a function of J3, 
for various fixed values of J2 in the range 0.4 < J3 < 1.0, 
for the anti-Neel and the spiral states of the spin-i J1-J2-J3 
model on the honeycomb lattice (with Ji = 1). The extrap- 
olated LSUBm (m — ► 00) results are based on the extrap- 
olation scheme of Eq. (|12|) . and the calculated results with 
m = {4, 6, 8} for both the aN and spiral phases for values 
J 2 = 0.69, 0.695, 0.7, 0.8, 0.9, 1.0; and with m = {6, 8, 10} for 
the aN phase for values J2 = 0.4, 0.5, 0.6. For the spiral state 
the results use the pitch angle <j> = ^lshb™ that minimizes the 
energy E — £?LSUBm (</>)• We note that in all cases curves with- 
out symbols attached refer to the anti-Neel state, whereas the 
corresponding curves with symbols refer to the spiral state. 



aN state to be shown only up to certain values of J3 for 
each curve shown in Fig. [TO] In each case, the LSUBm 
solution with a finite value of m extends further into the 
region where the aN solution actually ceases to exist (i.e., 
to after the energy crossing point with the spiral phase) . 
Presumably, in the m — > 00 limit, the LSUBto termi- 
nation points for the aN phase would coincide with the 
phase boundary with the spiral phase. Simple heuristic 
extrapolations based on the results with m = {4, 6, 8} 
agree well with the energy crossing points, and are hence 
entirely consistent with this hypothesis. 

We note from Fig. [10] that as the value of the pa- 
rameter Ji is decreased towards the lower limiting value 
J2 ~ 0.635, below which the spiral state ceases to exist, 
the energy curves for the aN and spiral phases lie increas- 
ingly close to one another, and hence the position of the 
crossing point becomes increasingly difficult to determine 
with high precision. Accordingly, we expect that a bet- 
ter indicator of the phase boundary might be obtained 
from a comparison of the magnetic order parameters of 
the two states, as now shown in Fig. [11] 

We see very clearly that for values of J2 ^ 0.69 the 
curves for the magnetic order parameters M of the aN 
and spiral phases cross in the physical regime (i.e., at a 
positive value of M). It is these crossing points that are 
our best estimates for the corresponding points on the 
boundary between the two phases, and these are shown 
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in the phase diagram of Fig. [2] by (magenta) times ( x ) 
symbols. For values J 2 > 0.8, these values are in excel- 
lent quantitative agreement with the corresponding en- 
ergy crossing points from Fig. [TDJ For smaller values of 
J2, down to the value J2 ~ 0.635 below which the spiral 
state ceases to exist for any value of J3, the energy cross- 
ing points become increasingly difficult to estimate accu- 
rately, as discussed above, and generally lie slightly below 
the much more accurate values obtained from Fig. [Til al- 
though they are still in good qualitative agreement with 
them. 

Our best estimate for the position of the third tricrit- 
ical QCP, (J 2 C3 , J3 3 ) = (0.69 ± 0.01,0.12 ± 0.02), which 
marks the point where the spiral and aN phases meet 
the QP phase, comes from curves such as those shown in 

Fig.rru 

For values J2 < Jj 3 ~ 0.69, we use the correspond- 
ing values of J3 at which the magnetic order parameter 
M — > for the aN phase, as shown in Fig. [Til to find the 
phase boundary between the aN and QP phases. The 
corresponding points are shown in the phase diagram of 
Fig. H] by (blue) plus (+) symbols. 



E. The quantum paramagnetic (PVBC?) phase(s) 

From the results presented so far we have seen that 
in the parameter space window J2, J3 E [0, 1] the spin-i 
Jy~J2~J?> HAFM on the honeycomb lattice with J\ = 1 
has regions of five different GS phases. Four of these 
(viz., the N, S, aN, and spiral phases) are quasiclassical 
in nature, and they almost completely surround the fifth 
QP phase, as shown in Fig. [21 with each of them sharing 
a boundary with the (almost) enclosed region of the QP 
phase. (Indeed, it seems likely that if the diagram were 
extended slightly to negative values of J3, the QP region 
would be seen to be entirely enclosed.) In the window 
J2tH G [0,1] there are three tricritical QCPs (and it 
seems likely that a fourth, which marks the meeting of 
the N, aN and QP phases, will occur just outside the 
window) . 

From our current results discussed here the question 
still remains open, however, as to the exact nature of 
this phase. What we know from previous wor k 15 ' 21 that 
employed the same CCM techniques as here is that the 
QP phase appears to be have PVBC ordering at least at 
four points along its boundary. These include the two 
points marked with the larger (red) times (x) symbols 
in Fig. [2] along its boundary with the N state where it 
crosses both the J3 = axis and the J3 = J2 line, the 
point marked with the larger (green) plus (+) symbol 
along its boundary with the S state where it crosses the 
J3 = J2 line, and the point marked with the larger (blue) 
plus (+) symbol along its boundary with the aN state 
where it crosses the J3 = axis. 

Those points were identified as lying on a phase bound- 
ary with the PVBC state by calculating, within the same 
CCM LSUBm approximations as used to calculate the 



phase transition points that marked the vanishing of the 
magnetic order parameter M in each case (i.e., for the N, 
S, and aN phases respectively), the susceptibility of the 
respective phases against the formation of PVBC order. 
We showed, within the accuracy of our results, that each 
of the above four points where the respective magnetic 
order parameter of each quasiclassical phase goes to zero 
coincide with the points at which the corresponding sus- 
ceptibility of the state to PVBC order becomes infinite. 

In principle we could now repeat those calculations for 
the PVBC susceptibility parameter for all points along 
the phase boundary of the QP state with the four qua- 
siclassical states. However, that would be particularly 
costly of computing resources for the spiral state. Even if 
we were to do so and hence show that the entire bound- 
ary region has PVBC order, we could still not be sure 
that the QP region was entirely PVBC-ordered, since 
there might still exists regions of other phases with other 
forms of order, possibly even of an exotic (spin-liquid) va- 
riety. The full characterization of the ordering within the 
QP phase(s) bounded by the four phases with quasiclas- 
sical ordering is an extremely challenging problem, and 
one that is essentially outside the scope of the present 
investigation. 



VI. SUMMARY AND DISCUSSION 

In this paper we have studied the spin-i HAFM on the 
honeycomb lattice with NN, NNN, and NNNN exchange 
interactions, namely the so-called J1-J2-J3 model de- 
scribed by the Hamiltonian of Eq. (TTJ). We have in- 
vestigated the full phase diagram of the model, in the 
case where all the bonds are antiferromagnctic in nature 
(i.e., J n > 0, n = 1,2,3), based on a combination of 
CCM techniques described in detail in Sees. IIIII and EI 
In particular we have set J\ = 1 to set the overall en- 
ergy scale, and we have restricted attention here to the 
window J2, J3 G [0, 1] for the remaining parameters. 

Our results are summarized in the phase diagram of 
Fig. [21 In the window J2, J3 G [0, 1] we find the five stable 
GS phases shown. Four of them are quasiclassical in na- 
ture, in the sense that they have counterparts in the clas- 
sical version (s — > 00) of the model. These comprise three 
phases showing collinear AFM order, viz., the Neel (N), 
striped (S), and anti-Neel states depicted in Figs. UJa), 
(b), and (d) respectively, plus a noncollinear spiral phase 
depicted in Fig. [TJc). In the classical version of the 
model, however, only the three phases with N, S, and 
spiral order exist in the examined window J2, J3 G [0, 1]. 
In the classical model the phase with aN order exists only 
for a part of the phase space where J3 < 0. We find that 
quantum fluctuations tend to stabilize this collinear aN 
phase at the expense of the spiral phase, in keeping with 
the very general observation that quantum fluctuations 
always seem to favor collinear phases over noncollinear 
ones. The fifth phase shown in Fig. [5] is a magnetically 
disordered, or quantum paramagnetic (QP), phase that 
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has no classical counterpart. 

The QP phase has phase boundaries with each of the 
four quasiclassical phases, with three tricritical QCPs oc- 
curring in the window J2, J3 G [0, 1], and a fourth pre- 
sumably occurring just outside the window with a small 
negative value of J3. The boundaries of the QP phase 
with each of the S, aN, and spiral phases appear to de- 
limit first-order transitions, whereas there are strong in- 
dications that the boundary of the QP phase with the 
N phase delimits a continuous phase transition. At two 
points along this boundary, viz., along the lines J3 = 
and J3 = J2; there arc strong indications that the QP 
phase there has PVBC order. Since the N and PVBC 
phases break different symmetries we have argued that 
the transitions there favor the dcconfincmcnt scenario. 
Nevertheless, we cannot entirely exclude even at these 
points the possibilities of a very weak first-order transi- 
tion or that the transition proceeds through a very nar- 
row region of intervening phase (possibly even of an ex- 
otic spin- liquid variety). 

We have found that all of the three remaining phase 
transition lines in the J2,Js G [0,1] window between 
the pairs of quasiclassical states shown in Fig. [5] (viz., 
between the N and S, S and spiral, and the spiral and 
aN phases) are first-order in nature. Rather strikingly, 
quantum effects in the s = | model turn the continuous 
transition between the spiral and S states of the classical 
(s — > 00) model into a first-order one. 

Whereas in the classical model the only three phases 
present in the J2,Js G [0,1] window (viz., the N, S, 
and spiral phases) meet at a single tricritical point at 
(J? ,J£' ) = (0.5,0.5), there are now three tricritical 
QCPs in the same window for the s = 5 model. The clas- 
sical tricritical point separates into two tricritical QCPs 
at (J 2 Cl , J 3 Cl ) = (0.51 ±0.01, 0.69 ±0.01) between the N, 
S, and QP phases, and at (J 2 2 , J 3 C2 ) = (0.65 ±0.02, 0.55 ± 
0.01) between the S, spiral, and QP phases. A third tri- 
critical QCP at (J 2 C3 , J 3 C3 ) = (0.69 ± 0.01,0.12 ± 0.02), is 
identified for the spin-i model between the spiral, aN, 
and QP phases. 

In overall terms our results for the phase diagram are 
in good agreement with other very recent studies of this 
model. For example, a study using a combination of vari- 
ous exact diagonalization (ED) and self-consistent cluster 
mean- field (SCCMF) techniques^ finds a phase diagram 
with basically the same five phases as we identify in the 
window J2, J3 G [0, 1]. There is good agreement with the 
positions of the two tricritical QCPs at (J^ 1 ,^ 1 ) an d 
(J 2 2 , J^ 2 )- The main difference seems to be in the posi- 
tion of the third tricritical point at (J 2 3 , J^ 3 )- Although 
both sets of calculations seem to be in good agreement 
about the boundary of the aN phase, the ED results gen- 
erally place the boundary between the spiral and QP 
phases to lower values of J2 such that the spiral phase 
occupies a larger region of phase space than our own 
calculations indicate. We note, however, that ED cal- 
culations are especially difficult for noncollincar phases, 
since the finite lattices used do not so easily sample such 



noncollinear phases. 

A further recent study of the model, using an un- 
biased pseudo-fcrmion functional renormalization group 
(PFFRG) method^ gives a phase diagram again in good 
overall agreement with ours, and with a phase boundary 
between the spiral and QP phases now in closer agree- 
ment with ours than from the ED results— Again, there 
is also good agreement with the phase boundaries in- 
volving the N and S states, and the positions of the two 
tricritical QCPs at (J 2 Cl , J 3 Cl ) and (J 2 2 , J 3 C2 ). The only 
qualitative disagreement is that the PFFRG study finds 
no evidence in the J2, J3 G [0, 1] window for the aN phase. 
Nevertheless, this study did find that for larger values of 
J2 and small values of J3 there were large incommen- 
surability shifts from the spiral phase and, furthermore, 
that there was evidence for J 2 > 0.4 of sizeable staggered 
dimer order. Such staggered dimer order is often difficult 
in practice to distinguish from the striped AFM order. 

Series expansion (SE) techniques have also been ap- 
plied to this model recently— In particular, expansions 
were performed around our N, S, and spiral states, as well 
as about the second classical spiral phase and the stag- 
gered dimer valence-bond crystal (SDVBC) state, which 
is also called the lattice nematic state. The results from 
the SE analysis are much more qualitative than ours or 
those of the ED+SCCMFi£ and PFFRG— studies. Nev- 
ertheless, the SE study is also in broad agreement, with 
the exception again of finding no evidence for the aN state 
(which might have been located for the only value J3 = 
that was studied by those authors with the second spiral 
phase that exists classically in the range g < J 2 < | for 
J3 = 0, and for which the aN phase is a limiting collincar 
form as discussed in Sec. ITT|) . 

It is clear that the existence or not of the aN phase in 
the J2, J3 G [0, 1] window is one of the points on which 
there is still disagreement between various studies. It 
seems clear that, if it does indeed exist, as we argue here, 
it becomes unstable at very small values of J3 for all 
values of J2 < 1 for which it exists. We note that a 
quite different recent ED study-ii of the model (along 
the J3 = line only), which used a NN singlet valence- 
bond basis, found very similar critical points to ours for 
the boundaries of the QP phase, but found that while the 
QP phase was bounded on one side (for smaller values of 
J 2 ) by the N phase, it was bounded on the other side by 
a SDVBC phase. However, these results would again be 
equally consistent with the identification of this phase as 
our aN phase, since their interpretation is strongly biased 
by their choice of basis. 

A quite separate study of the case J3 — has also 
been performed recently based on an entangled-plaqucttc 
variational (EPV) ansatz— This EPV study uses a sin- 
gle very broad class of entangled-plaquette states as trial 
wave functions, and finds in this very unbiased way that 
along the J3 = line the model has N, QP, and aN phases 
with critical points very close to ours. In the same re- 
gion studied by us (viz., J2 < 1) the EPV study finds no 
evidence of spiral order along the J3 = line. 
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It is clear that the spin-i J1-J2-J3 HAFM on the hon- 
eycomb lattice is a challenging model, but one in which 
there seems now to be a growing consensus on its over- 
all phase structure. There is very good agreement over 
the regions in which the N and S phases exist, and we 
believe our own CCM results now give perhaps the best 
quantitative results in these cases for the positions of the 
phase boundaries and the positions of the two tricritical 
QCPs at (J 2 C1 , J 3 C1 ) and (J| 2 , J 3 C2 ). 

An uncertainty remains over the precise extent of the 
phase with spiral order, and whether or not there is an 
aN phase along the J3 = line for values of J2 beyond 
the point where the QP phase disappears, and hence also 
for small positive values of J3 up to the point where spiral 
order sets in. The other main uncertainty is the nature 
of the QP phase itself. We have argued here that over 
at least some widely separated points on the boundaries 
with the N, S, and aN phases the QP phase has PVBC 
order. Two quite separate ED calculation a 16 ' 17 also give 
clear evidence that much of the QP phase has PVBC 
order, although the latter calculations^ are only done 
along the J3 = line. 

By contrast the EPV calculations^ along the J3 = 
line seem to favor a disordered (spin-liquid) phase, while 
spin-wave calculations^ favor SDVBC order along the 
same line in the QP regime. The PFFRG study^ also 
done over the entire J2, J3 € [0, 1] window, finds evidence 
too that a large part of the QP regime has strong SDVBC 
order, while the part with smaller values of J2 has only 
weak PVBC order. On the other hand the SE studyi£ 
finds that SDVBC order is not favored, at least for low 



values of J3. 

Wc should note, however, that the SE study is in broad 
disagreement with most other studies along the J3 = 
line, in that it finds no evidence at all for a magnetically 
disordered phase there, but instead finds that the N phase 
first gives way to the second classical spiral phase, and 
then later to the spiral phase considered here, as J 2 is in- 
creased. On the other hand, the SE results are consistent 
with the finding from the ED+SCCMF analysis^ that at 
least for some parameter ranges the SDVBC state might 
be very difficult to distinguish from magnetically ordered 
states such as our S state. 

Finally wc note that the ED+SCCMF study also 
presents evidence for the N to PVBC transition being 
a strong candidate for a deconfined transition, just as 
we have found in our earlier CCM studies of the model 
along the J3 = J2 linei£ and along the J 3 = line. 21 
Clearly this model still has open questions, but we be- 
lieve that the CCM results presented here have furthered 
our understanding of it. 
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